% This is a test code for the unified pole generation algorithm for
% Fermi-Dirac function. The Hamiltonian is taken to be a 2D Laplacian plus a
% small random external potential.
% 
% Inputvars:
% Npole: The number of poles to be calculated.
% beta: inverse temperature. In this code the unit is inverse hartree.
% deltaE: the spectrum width
% deltagap: Energy gap defined by min(abs(EV-mu))



Npole = 40;
beta = 1052;
deltaE = 5;
deltagap = 0.00001;


display('Generating Hamiltonian...');
Npolehalf = Npole/2;
if(1)
	N = 16;
	h = 1/N;
	epsilon = 1 / N;
	effh = sqrt(2) * h / epsilon; 
	% effh = h / epsilon; 
	% V = (1d-1*randn(N/8,N/8) - 0.2391);
	% V = interpft(interpft(V,N,1),N,2);

	NSite = (N/8)^2+1;
	NElectron = NSite * 2;
	mu = 0.0953;
	if( 1 ) 
		V = 0.001*rand(N, N);
		save ranpot V;
	else
		load ranpot V;
	end

	isRhoComput = 1;
	isEgyComput = 0;
	isMuUsed = 0;
	[H, RhoExact, EgyExact, mu1] = ham2d5( N, V, mu, effh, beta, ...
		isRhoComput, isEgyComput, isMuUsed, NSite );
end

if( isMuUsed == 0 )
	mu = mu1;
end

disp('Calculating poles...')

M = deltaE;
mshift = (pi/beta)^2;
m2 = mshift+(deltagap)^2;
M2 = M^2;
k = (sqrt(M2/m2)-1)/(sqrt(M2/m2)+1);
L = -log(k)/pi;
[K,Kp] = ellipkkp(L);


t = .5i*Kp - K + (.5:Npolehalf)*2*K/Npolehalf;  %t = .9999999999i*Kp - K + (.5:Npolehalf)*2*K/Npolehalf;
[u cn dn] = ellipjc(t,L);
z = sqrt(m2*M2)*((1/k+u)./(1/k-u));
dzdt = cn.*dn./(1/k-u).^2;

zsqrt = sqrt(z-mshift);

zall = [zsqrt -zsqrt];

zweight = zeros(Npole, 1);
zshift = zeros(Npole, 1);

for j = 1 : Npolehalf
	zshift(j) = -mu - zsqrt(j);
	zshift(j+Npolehalf) = -mu + zsqrt(j);
	zweight(j) = tanh(beta*zsqrt(j)/2)*(-1) / zsqrt(j) * dzdt(j) * ...
		2*K*sqrt(m2*M2)/(k*pi*Npolehalf);
	zweight(j+Npolehalf) = zweight(j);
end

display('Calculating...');

S = ones(N, N);
for j = 1:length(zweight)
	j
	Vnow = V + zshift(j);
	Res = five_extractall(Vnow, N, effh);
	S = S + imag(zweight(j)*Res);
end       
errrho = sum(sum(abs(S-RhoExact))) / NElectron


display('Outputing poles...');

fid = fopen( 'pole.in', 'w' );
fprintf( fid, '%5i\n', Npole ) ;
for j = 1 : length(zshift)
	fprintf( fid, '(%25.15e,%25.15e)\n', real(zshift(j)), imag(zshift(j)) ); 
end
for j = 1 : length(zweight)
	fprintf( fid, '(%25.15e,%25.15e)\n', real(zweight(j)), imag(zweight(j)) );
end
fclose(fid);
